
Job: betamc.std.prg (/home/gov/faculty/ppaolino/parr)
Sun Jun 17 10:33:17 2001, pid=899, uid=1141(ppaolino)

cheng alg used if cheng=1, johnk otherwise, cheng=       0.0000000 
parm1 =      -0.80000000       0.30000000 
parm2 =      -0.40000000      -0.30000000 
p =       0.46193747 
q =       0.72968919 
true first difference (mean)=      0.31021425 
true first difference (variance)=     0.036388892 

descriptive stats for standard beta parameters

     -0.76548828       0.11629184 
      0.30718160       0.11785654 
     -0.36244088       0.13950947 
     -0.30140043       0.13036951 
 
standard beta mse=     0.099048853 
 
alternative beta mse=     0.099189029 
 
ols mse=     0.099097400 
 
logit trans ols mse=      0.12030862 
 
het ols mse=     0.099171002 
 
logit trans het ols mse=      0.11870367 
 
log trans het ols mse=      0.13054685 
 
first differences for mean and variance for std beta (sim)
      0.30924087      0.061422933 
     0.033389675      0.017661418 
 
first differences for mean and variance for std beta (nosim)
      0.31168599      0.061802352 
     0.034656035      0.018057309 
 
first differences for mean and variance for alt beta (sim)
      0.29056825      0.054970888 
     0.031701608      0.017541176 
 
first differences for mean and variance for alt beta (nosim)
      0.29273687      0.055251639 
     0.032780519      0.017843375 
 
first differences for mean ols
      0.29358494      0.057546501 
 
first differences for mean ols (logit trans)
      0.56799893       0.11446788 
 
first differences for mean and variance for het normal
      0.28782282      0.054306853 
     0.035425857      0.026935100 
 
first differences for mean and variance for logit trans.
      0.54929139       0.11375074 
      0.37699428       0.15327127 
 
first differences for mean ols (log trans)
      0.26912355      0.063909688 
 
efficiency of standard beta vs. alternative beta=       93.744322 
 
efficiency of beta vs. ols=       96.899788 
 
efficiency of beta vs. ols (log trans)=       456.44587 
 
efficiency of standard beta vs. het normal=       95.027939 
 
efficiency of standard beta vs. het normal (log trans)=       428.44886 
 
efficiency of standard beta vs. OLS (log trans)=       122.92270 
 
efficiency of standard beta vs. alternative beta=       100.35601 
 
efficiency of standard beta vs. het normal=       148.57669 
 
efficiency of standard beta vs. het normal (log trans)=       2059.8168 
 
standard beta parms' overconfidence=
       92.515376 
       102.78019 
       100.18681 
       100.08573 
 
alternative beta parms' overconfidence=
       98.723983 
       97.032287 
       97.052287 
       97.497043 
 
ols parms' overconfidence=
       96.955268 
       89.985218 
 
ols (log. trans) parms' overconfidence=
       91.833117 
       122.30827 
 
ols (log trans) parms' overconfidence=
       89.135602 
       94.070956 
 
heteroskedastic normal parms' overconfidence=
       98.148485 
       85.999960 
       71.762877 
       85.227630 
 
heteroskedastic log. trans. parms' overconfidence=
       93.320420 
       117.35533 
       149.81695 
       166.96971 
 
standard beta's mean first difference overconfidence
       102.70816 
 
standard beta's variance first difference overconfidence
       98.263925 
 
alternative beta's mean first difference overconfidence
       97.255456 
 
alternative beta's variance first difference overconfidence
       96.579576 
 
check standard beta's percentage within 95% conf int

      0.96200000 
      0.94600000 
      0.95000000 
      0.94800000 
 
check alternative beta's percentage within 95% conf int

      0.95800000 
      0.96300000 
      0.95500000 
      0.95300000 
 
check ols percentage within 95% conf int

      0.95600000 
      0.97100000 
 
check logit trans/ols percentage within 95% conf int

      0.96400000 
      0.88500000 
 
check het ols percentage within 95% conf int

      0.95400000 
      0.97500000 
      0.99200000 
      0.98200000 
 
check logit trans/ols percentage within 95% conf int

      0.96200000 
      0.91000000 
      0.80500000 
      0.75600000 
 
check logit trans/ols percentage within 95% conf int

      0.96800000 
      0.96600000 
 
check percentage of sig parms (standard beta)
      0.99700000 
      0.47400000 
 
check percentage of sig parms (alternate beta)

      0.87000000 
      0.99700000 
      0.33000000 
     0.075000000 
 
check percentage of sig parms (basic ols)

       1.0000000 
      0.99600000 
 
check percentage of sig parms (logit trans ols)

      0.85300000 
      0.99700000 
 
check percentage of sig parms (het ols)

       1.0000000 
      0.99600000 
       1.0000000 
      0.14800000 
 
check percentage of sig parms (het ols - logit trans)

      0.86500000 
      0.99500000 
       1.0000000 
      0.38000000 
 
check percentage of sig parms (ols -- log tranform)

       1.0000000 
      0.99500000 
 
descriptive stats for y's for successful runs

      0.39891303 
      0.34970059 
   2.1624152e-05 
      0.99881989 
descriptive stats for y's for failed runs
       0.0000000 
descriptive stats for x's for successful runs

    -0.093291739 
       1.1236755 
      -3.3593560 
       2.4543428 
descriptive stats for x's for failed runs
       0.0000000 
number of failed runs
       0.0000000 

code used to generate the above results:

library cml;
#include cml.ext;
CMLset;

fails=0;
crash=0;
ystatf=0;
xstatf=0;

/* load x matrix based upon number of observations */

x = ones(100,1)~rndn(100,1);
n = rows(x);

/* produce matrices for simulated first differences */

meanx=meanc(x); @ create vector of means for the indep vars @
stdx=stdc(x);   @ create vector of stdev for the indep vars @
minx=minc(x);
maxx=maxc(x);

/* start obtaining the first difference vectors (means) */

sw=zeros(1,cols(x)-1)|eye(cols(x)-1);
xl=meanx-stdx.*sw;           @ these are kXk-1 vectors @
xh=meanx+stdx.*sw;       
xmin=meanx.*(1-sw)+minx.*sw;  
xmax=meanx.*(1-sw)+maxx.*sw;

/* start creating Monte Carlo data */
k = cols(x);
parm1 = {-0.8, 0.3};
parm2 = {-0.4, -0.3};
cheng=0;
print "cheng alg used if cheng=1, johnk otherwise, cheng=";;cheng; 
p = exp(x*parm1);
q = exp(x*parm2);

/* determine true first differences */

/* get  alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         alt=exp(xl'*parm1);                          @ nxk-1 matrix @
         blt=exp(xl'*parm2);
         aht=exp(xh'*parm1);
         bht=exp(xh'*parm2);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlt=alt./(alt+blt);                    @ nxk-1 matrix @
         meanht=aht./(aht+bht);
         varlt=(alt.*blt)./(((alt+blt)^2).*(alt+blt+1));
         varht=(aht.*bht)./(((aht+bht)^2).*(aht+bht+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmt=meanht-meanlt;
         fdvt=varht-varlt;

/* print data parameters for simulation */

print "parm1 = ";; parm1';
print "parm2 = ";; parm2';
print "p = ";; meanc(exp(x*parm1));
print "q = ";; meanc(exp(x*parm2));
print "true first difference (mean)=";; fdmt;
print "true first difference (variance)=";; fdvt;

/* start Monte Carlo trials */
trial=1;
do while trial<=1000; 

print "MC trial=";; trial;

   if trial==1;
      ret=0;
      if fails==0;
         bt=0; 
      endif;
   endif;

   if fails==0;
   success=0;

/* comment out random number generator if p,q<1 */
if cheng==1;
/* create J&K (Chen) alpha */
ralpha=p+q;

/* create J&K (Chen) beta */
p1=p.>1;
q1=q.>1;
rbeta=real((p1.*q1).*(sqrt((ralpha-2)./(2.*p.*q-ralpha))) +
      (1-p1.*q1).*(1/(minc((p'|q')))));

/* create J&K (Chen) gamma */
rgamma=p+1/rbeta;

/* step 1 in J&K (Chen) RNG */
u1=rndu(n,1);
u2=rndu(n,1);
v=rbeta.*ln(u1./(1-u1));
w=p.*exp(v);

/* step 2 in J&K (Chen) RNG */

done=0;
      do while done<n;
         d=(ralpha.*ln(ralpha./(q+w))+rgamma.*v-1.3862944).<ln((u1^2).*u2); 
         u1=(1-d).*u1 + d.*rndu(n,1);
         u2=(1-d).*u2 + d.*rndu(n,1);
         v=(1-d).*v + d.*rbeta.*ln(u1./(1-u1));
         w=(1-d).*w + d.*p.*exp(v);
         done=n-sumc(d);
      endo;

/* step 3 in J&K (Chen) RNG */
y=w./(q+w); 
endif;

/* discretize y 

y1=yc.>0;
y2=yc.>0.14;
y3=yc.>0.28;
y4=yc.>0.42;
y5=yc.>0.56;
y6=yc.>0.70;
y7=yc.>0.84;

y=(y1+y2+y3+y4+y5+y6+y7)/8; */

/* comment out random number generator if p or q>1 */
if cheng==0;
/* step 1 in Johnk's RNG */

u1=rndu(n,1);
u2=rndu(n,1);

v1=u1^(1/p);
v2=u2^(1/q);

/* step 2 in Johnk's RNG */

w=v1+v2;

/* step 3 in Johnk's RNG */

done=0;
      do while done<n;
         d=w.>1;
         u1=(1-d).*u1 + d.*rndu(n,1);
         u2=(1-d).*u2 + d.*rndu(n,1);
         v1=(1-d).*v1 + d.*u1^(1/p);
         v2=(1-d).*v2 + d.*u2^(1/q);
         w=v1+v2;
         done=n-sumc(d);
      endo;

/* step 4 in Johnk's RNG */
y=v1./w;
endif;

/* cheat function to re-set Ys that equal 1 --- would choke llik function
   otherwise */
d=y.==1;
print "cheat=";;sumc(d);
y=(1-d).*y + d.*.99999;  


/* return the monte carlo data set */
    data = y~x[.,2:cols(x)];

   endif;

   if fails==0;

      do while success==0 and fails<=30;

/* maximum likelihood function for estimating parameters alpha and beta */
__title=" ** A Beta Model, Analytic Solution (Standard Estimation) **";

stval=1.0|ones(k-1,1)*0.5|1.0|ones(k-1,1)*0.4;

/* for failed runs, use starting values of crashed iteration */

       if fails>=1 and fails<=40;
          stval=(bt')/2;
          print "retrying run";
       endif;

print "mean y, max y, min y, std y";
meanc(y)~maxc(y)~minc(y)~stdc(y);

proc psi(x);
  local p;
  x=x+6;
  p=1/(x.*x);
  p=(((0.004166666666667*p-0.003968253986254).*p+
      0.008333333333333).*p-0.083333333333333).*p;
  p=p+ln(x)-0.5/x-1/(x-1)-1/(x-2)-1/(x-3)-1/(x-4)-1/(x-5)-1/(x-6);
retp(p);
endp;

proc gradproc(theta,ds);
    local y,x,z,be,h,xb,p,q,per1,gr1,gr2,grad;
		     
    y=ds[.,1];
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    be=theta[1:cols(x),1];                 @ parameters for alpha @
    h=trimr(theta,rows(be),0);            @ parameters for beta @
    p=exp(x*be);                 @ expression for alpha @
    q=exp(x*h);                  @ expression for beta @

    per1=p+q;
    gr1= psi(p+q).*x.*p-psi(p).*x.*p+x.*p.*ln(y);
    gr2= psi(p+q).*x.*q-psi(q).*x.*q+x.*q.*ln(1-y);
    grad=gr1~gr2;
    retp(grad);
    endp;

proc loglik(theta,ds);
    local y,x,z,be,h,xb,p,q,l1,l2,l3,l4,l5,llik;
		     
    y=ds[.,1];
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    be=theta[1:cols(x),1];                 @ parameters for alpha @
    h=trimr(theta,rows(be),0);            @ parameters for beta @
    p=exp(x*be);                 @ expression for alpha @
    q=exp(x*h);                  @ expression for beta @

    if (sumc((p.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((q.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((y.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc(((1-y).>0)))<n;
       llik=miss(0,0); print "run failed"; 
    else;
       llik=lngm(p+q)-lngm(p)-lngm(q)+(p-1).*ln(y)+
         (q-1).*ln(1-y);  @ log-likelihood @
    endif;
    retp(llik);endp;

_cml_GradProc=&gradproc;
_cml_GradCheck=1;
_cml_Options = { none }; 
_cml_MaxIters=400;

   if fails.>=5 and fails.<=10;
     _cml_DirTol=.000001;
   endif;
   if fails.>10 and fails.<=15;
     _cml_DirTol=.0000001;
   endif;
   if fails.>15 and fails.<=20;
     _cml_DirTol=.00000001;
   endif;
   if fails.>20;
     _cml_DirTol=.00001;
   endif;

    {b,logl,g,vc,ret}=CMLprt(CML(data,0,&loglik,stval));

/* check that parameters are reasonable */
   bt=(b');

/* determine success or failure of the trial */
   cor=vc[1,1]/vc[1,1];
         if ret>=1 or cor/=1;
            fails=fails+1;
         elseif ret==0 and cor==1;
            success=1; fails=0;
         endif;
      endo;
   endif;

   if fails==0;
      success=0; bta=0;
      do while success==0 and fails<=30;

CMLset;
/* Alternative beta estimation */
__title=" ** A Beta Model, Numerical Solution (Alternative Estimation) **";

stvala=0.0|ones(k-1,1)|2.0|ones(k-1,1)*0.4;

/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=30;
         stvala=(bta')/2;
         print "retrying run";
      endif;

proc gradpa(theta,ds);
    local y,one,x,z,be,h,xb,zh,zh1,e,per1,gr1,gr2,grad;
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=exp(x*be);
    zh=exp(z*h);
    zh1=zh+1;
    e=xb./(1+xb);
    per1=(-x.*e+x.*(e^2)).*zh1;

    gr1=psi(e.*zh1+(1-e).*zh1-1).*(x.*e.*zh1-(e^2).*zh1.*x+per1)-
        psi(e.*zh1-e).*(x.*e.*zh1-(e^2).*zh1.*x-x.*e+(e^2).*x)-
        psi((1-e).*zh1-1+e).*(per1+x.*e-(e^2).*x)+
        (x.*e.*zh1-(e^2).*zh1.*x-x.*e+(e^2).*x).*ln(y)+
        (per1+x.*e-(e^2).*x).*ln(1-y);
    gr2=psi(e.*zh1+(1-e).*zh1-1).*(e.*z.*zh+(1-e).*z.*zh)-
        (psi(e.*zh1-e).*xb.*z.*zh)./(1+xb)-
        psi((1-e).*zh1-1+e).*(1-e).*z.*zh+xb.*z.*zh.*ln(y)./(1+xb)+
        (1-e).*z.*zh.*ln(1-y);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglika(theta,ds);
    local y,one,x,z,be,h,xb,zh,e,v,p,q,llik,ds1;
		     
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=exp(x*be);
    zh=exp(z*h);
    e=xb./(1+xb);
    v=e.*(1-e).*(1/(zh+1));
    p=((e^2).*(1-e))./v-e;
    q=(e.*(1-e)^2)./v-(1-e);
    if (sumc((p.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((q.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((y.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc(((1-y).>0)))<n;
       llik=miss(0,0); print "run failed"; 
    else;
       llik=lngm(p+q)-lngm(p)-lngm(q)+(p-1).*ln(y)+
         (q-1).*ln(1-y);  @ log-likelihood @
    endif;
    retp(llik);endp;

_cml_GradProc=&gradpa;
_cml_GradCheck=1;
_cml_Options = { none }; 
_cml_MaxIters=400; 
    {ba,logl,g,vca,ret}=CMLprt(CML(data,0,&loglika,stvala));

/* check that parameters are reasonable */
   bta=(ba');

/* determine success or failure of the trial */
   cora=vca[1,1]/vca[1,1];
         if ret>=1 or cora/=1;
            fails=fails+1;
         elseif ret==0 and cora==1;
            success=1; fails=0;
         endif;
      endo;
   endif;

/* linear heteroskedastic estimation */

   if fails==0;
      success=0; bth=0;
      do while success==0 and fails<=30;

CMLset;

__title=" ** A Linear-Normal Model, Numerical Solution **";

stvalh=inv(x'*x)*x'*y|(-0.5*ones(cols(x),1));
       
/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=10;
         stvalh=(bth')/2;
         print "retrying run";
      endif;  

proc gradph(theta,ds);
    local y,one,x,z,be,h,xb,zh,gr1,gr2,grad;
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=x*be;
    zh=exp(z*h);
    gr1=(y-xb).*x./zh;
    gr2=-0.5*(x-(x.*(y-xb)^2)./zh);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglikh(b,ds);                     
    local xb,x1,y,llik,b1,g1,s,s2,u,ds1,one;
    y=ds[.,1];    /* dependent variable */
    one=ones(rows(ds),1);
    x1=one~ds[.,2:cols(ds)];
    b1=b[1:cols(x1),1];   /* partition off parms for E(y) */
    xb=x1*b1;
    s=trimr(b,rows(b1),0);      /* partition off parms for V(e) */
    g1=one~ds[.,2:cols(ds)];
    s2=exp(g1*s);
    u=y-xb;             
    llik=-0.5*(ln(s2)+u.*u./s2);  /* likelihood function */
    retp(llik);endp;

_cml_GradProc=&gradph;
_cml_GradCheck=1;
_cml_Options = { none };
_cml_DirTol=0.000001; 
_cml_MaxIters=500;

    {blh,logl,g,vch,ret}=CMLprt(CML(data,0,&loglikh,stvalh));

/* check that Monte carlo parameters are reasonable */
bth=(blh');

/* determine success or failure of the trial */
   corh=vch[1,1]/vch[1,1];
       if ret>=1 or corh/=1;
          fails=fails+1;
          elseif ret==0 and corh==1;
          success=1; fails=0;
       endif;
    endo;
endif;

/* linear heteroskedastic estimation (with logit tranformation) */

   if fails==0;
      success=0; btl=0;
      do while success==0 and fails<=30;

CMLset;

__title=" ** A Linear-Normal Model, Numerical Solution **";

stvall=inv(x'*x)*x'*(ln(y./(1-y)))|(-0.5*ones(cols(x),1));
       
/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=30;
         stvall=(btl')/2;
         print "retrying run";
      endif;  

proc gradpl(theta,ds);
    local y,one,x,z,be,h,xb,zh,gr1,gr2,grad;
    y=ln(ds[.,1]./(1-ds[.,1]));
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=x*be;
    zh=exp(z*h);
    gr1=(y-xb).*x./zh;
    gr2=-0.5*(x-(x.*(y-xb)^2)./zh);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglikl(b,ds);                     
    local xb,x1,y,llik,b1,g1,s,s2,u,ds1,one;
    y=ln(ds[.,1]./(1-ds[.,1]));    /* dependent variable */
    one=ones(rows(ds),1);
    x1=one~ds[.,2:cols(ds)];
    b1=b[1:cols(x1),1];   /* partition off parms for E(y) */
    xb=x1*b1;
    s=trimr(b,rows(b1),0);      /* partition off parms for V(e) */
    g1=one~ds[.,2:cols(ds)];
    s2=exp(g1*s);
    u=y-xb;             
    llik=-0.5*(ln(s2)+u.*u./s2);  /* likelihood function */
    retp(llik);endp;

_cml_GradProc=&gradpl;
_cml_GradCheck=1;
_cml_Options = { none };
_cml_DirTol=0.000001; 
_cml_MaxIters=500;

    {bll,logl,g,vcl,ret}=CMLprt(CML(data,0,&loglikl,stvall));

/* check that Monte carlo parameters are reasonable */
btl=(bll');

/* determine success or failure of the trial */
   corl=vcl[1,1]/vcl[1,1];
      if ret>=1 or corl/=1;
         fails=fails+1;
      elseif ret==0 and corl==1;
         success=1; fails=0;
      endif;
   endo;
endif;

/* if trial has failed 30 times */

   if fails==30;
      print "run failed after 30 attempts";
      crash=crash+1;

      if ystatf==0;
         ystatf=(meanc(y)~stdc(y)~minc(y)~maxc(y));
      elseif ystatf/=0;
         ystatf=ystatf|(meanc(y)~stdc(y)~minc(y)~maxc(y));
      endif;

      if xstatf==0;
         xstatf=(meanc(x[.,2])|stdc(x[.,2])|minc(x[.,2])|maxc(x[.,2]));
      elseif xstatf/=0;
         xstatf=xstatf~(meanc(x[.,2])|stdc(x[.,2])|minc(x[.,2])|
                 maxc(x[.,2]));
      endif;
      fails=0;
      clear y;
   endif;

/* collect data for a successful trial */
   if success==1;

/* collect y stats for a successful trial */
      if trial==1;
         ystats=meanc(y)~stdc(y)~minc(y)~maxc(y);
         xstats=(meanc(x[.,2:cols(x)])|stdc(x[.,2:cols(x)])|
                minc(x[.,2:cols(x)])|maxc(x[.,2:cols(x)]));
         bsum=b;
         sesum=sqrt(diag(vc));
         bsuma=ba;
         sesuma=sqrt(diag(vca));
         bsumh=blh;
         sesumh=sqrt(diag(vch));
         bsuml=bll;
         sesuml=sqrt(diag(vcl));
      endif;
      if trial>=2;
         ystats=ystats|(meanc(y)~stdc(y)~minc(y)~maxc(y));
         xstats=xstats~(meanc(x[.,2:cols(x)])|stdc(x[.,2:cols(x)])|
                minc(x[.,2:cols(x)])|maxc(x[.,2:cols(x)]));
         bsum=bsum~b;
         sesum=sesum~(sqrt(diag(vc)));
         bsuma=bsuma~ba;
         sesuma=sesuma~sqrt(diag(vca));
         bsumh=bsumh~blh;
         sesumh=sesumh~sqrt(diag(vch));
         bsuml=bsuml~bll;
         sesuml=sesuml~sqrt(diag(vcl));
      endif;

/* vectors of first difference values for the j^th independent var */
         bns=b[1:cols(x),.];                @ nxk matrix @
         hns=b[cols(x)+1:rows(b),.];   @ nxk matrix @

/* calculate model mse */

         ahat=exp(x*bns);
	 bhat=exp(x*hns);
         yhatns=ahat./(ahat+bhat);
	 bnsmse=sumc((y-yhatns)^2)/n;
	 
/* get predicted alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         alns=exp(xl'*bns);                          @ nxk-1 matrix @
         blns=exp(xl'*hns);
         ahns=exp(xh'*bns);
         bhns=exp(xh'*hns);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlns=alns./(alns+blns);                    @ nxk-1 matrix @
         meanhns=ahns./(ahns+bhns);
         varlns=(alns.*blns)./(((alns+blns)^2).*(alns+blns+1));
         varhns=(ahns.*bhns)./(((ahns+bhns)^2).*(ahns+bhns+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmns=meanhns-meanlns;
         fdvns=varhns-varlns;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totmns=fdmns;  @ k-1X1 matrix for mean @
            totvns=fdvns;  @ k-1X1 matrix for variance @
	    totbmse=bnsmse;
         endif;

         if trial>=2 and ret==0; 
            totmns=totmns~fdmns;  @ k-1X(trial) matrix for mean @
            totvns=totvns~fdvns;  @ k-1X(trial) matrix for variance @
	    totbmse=totbmse|bnsmse;
         endif;

/* vectors of first difference values for the j^th independent var */
         betasim=(rndmn(b,vc,1000));  @ nxk matrix, where n=# of draws @
         be=betasim[.,1:cols(x)];                @ nxk matrix @
         h=betasim[.,cols(x)+1:cols(betasim)];   @ nxk matrix @

/* get predicted alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         al=exp(be*xl);                          @ nxk-1 matrix @
         bl=exp(h*xl);
         ah=exp(be*xh);
         bh=exp(h*xh);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlow=al./(al+bl);                    @ nxk-1 matrix @
         meanhigh=ah./(ah+bh);
         varlow=(al.*bl)./(((al+bl)^2).*(al+bl+1));
         varhigh=(ah.*bh)./(((ah+bh)^2).*(ah+bh+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdm=meanhigh-meanlow;
         fdv=varhigh-varlow;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totm=meanc(fdm);  @ k-1X1 matrix for mean @
            totv=meanc(fdv);  @ k-1X1 matrix for variance @
            totms=stdc(fdm);  @ k-1X1 matrix for mean @
            totvs=stdc(fdv);  @ k-1X1 matrix for variance @
         endif;

         if trial>=2 and ret==0; 
            totm=totm~meanc(fdm);  @ k-1X(trial) matrix for mean @
            totv=totv~meanc(fdv);  @ k-1X(trial) matrix for variance @
            totms=totms~stdc(fdm);  @ k-1X(trial) matrix for mean @
            totvs=totvs~stdc(fdv);  @ k-1X(trial) matrix for variance @
         endif;

/* get relevant stats for alternative method of beta estimation */
/* vectors of first difference values for the j^th independent var */
         bnsa=ba[1:cols(x),.];                @ nxk matrix @
         hnsa=ba[cols(x)+1:rows(ba),.];   @ nxk matrix @

/* calculate mse for alternative beta estimation */
  
         yhatnsa=exp(x*bnsa)./(1+exp(x*bnsa));
 	 bnsamse=sumc((y-yhatnsa)^2)/n;
	 
/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlnsa=(exp(xl'*bnsa)./(1+exp(xl'*bnsa)))';   @ nxk-1 matrix @
         meanhnsa=(exp(xh'*bnsa)./(1+exp(xh'*bnsa)))';
         dlsa=(1/(exp(xl'*hnsa)+1))';
         dhsa=(1/(exp(xh'*hnsa)+1))';
         varlnsa=meanlnsa.*(1-meanlnsa).*dlsa;
         varhnsa=meanhnsa.*(1-meanhnsa).*dhsa;

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmnsa=meanhnsa-meanlnsa;
         fdvnsa=varhnsa-varlnsa;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totmnsa=fdmnsa;  @ k-1X1 matrix for mean @
            totvnsa=fdvnsa;  @ k-1X1 matrix for variance @
	    totbamse=bnsamse;
         endif;

         if trial>=2 and ret==0; 
            totmnsa=totmnsa~fdmnsa;  @ k-1X(trial) matrix for mean @
            totvnsa=totvnsa~fdvnsa;  @ k-1X(trial) matrix for var @
	    totbamse=totbamse|bnsamse;
         endif;

/* vectors of first difference values for the j^th independent var */
         betasima=(rndmn(ba,vca,1000));  @ nxk matrix, where n=# of draws @
         bea=betasima[.,1:cols(x)];                @ nxk matrix @
         ha=betasima[.,cols(x)+1:cols(betasima)];   @ nxk matrix @

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanla=(exp(bea*xl)./(1+exp(bea*xl)));   @ nxk-1 matrix @
         meanha=(exp(bea*xh)./(1+exp(bea*xh)));
         dl=(1/(exp(ha*xl)+1));
         dh=(1/(exp(ha*xh)+1));
         varla=meanla.*(1-meanla).*dl;
         varha=meanha.*(1-meanha).*dh;

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdma=meanha-meanla;
         fdva=varha-varla;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totma=meanc(fdma);  @ k-1X1 matrix for mean @
            totva=meanc(fdva);  @ k-1X1 matrix for variance @
            totmsa=stdc(fdma);
            totvsa=stdc(fdva);
         endif;

         if trial>=2 and ret==0; 
            totma=totma~meanc(fdma);  @ k-1X(trial) matrix for mean @
            totva=totva~meanc(fdva);  @ k-1X(trial) matrix for variance @
            totmsa=totmsa~stdc(fdma);
            totvsa=totvsa~stdc(fdva);
         endif;

/* set parameter vectors (linear heteroskeadastic) */
         bhh=blh[1:cols(x),1];
         hh=trimr(blh,rows(bhh),0);

/* calculate mse */

         yhatlh=x*bhh;
	 lhmse=sumc((y-yhatlh)^2)/n;
	 
/* get predicted means */
         elh=xl'*bhh;
         ehh=xh'*bhh;

/* get predicted variances */
         vlh=exp(xl'*hh);
         vhh=exp(xh'*hh);

/* produce first differences */
         if trial==1 and ret==0;
         totfdmh=ehh-elh;
         totfdvh=vhh-vlh;
	 totlhmse=lhmse;
         elseif trial>=2 and ret==0;
         totfdmh=totfdmh~(ehh-elh);
         totfdvh=totfdvh~(vhh-vlh);
	 totlhmse=totlhmse|lhmse;
         endif;

/* set parameter vectors (linear logit tranformation) */
         blo=bll[1:cols(x),1];
         hl=trimr(bll,rows(blo),0);
	 
/* calculate mse for linear het log trans */
  
         yhatlhlt=exp(x*blo)./(1+exp(x*blo));
	 lhltmse=sumc((y-yhatlhlt)^2)/n;

/* get predicted means */
         eol=exp(xl'*blo)./(1+exp(xl'*blo));
         eoh=exp(xh'*blo)./(1+exp(xh'*blo));

/* get predicted variances */
         vol=((exp(xl'*blo)./(1+exp(xl'*blo))^2)^2).*exp(xl'*hl);
         voh=((exp(xh'*blo)./(1+exp(xh'*blo))^2)^2).*exp(xh'*hl);

/* produce first differences */
         fdml=eoh-eol;
         fdvl=voh-vol;

/* produce first differences */
         if trial==1 and ret==0;
         totfdml=fdml;
         totfdvl=fdvl;
	 totthmse=lhltmse;
         elseif trial>=2 and ret==0;
         totfdml=totfdml~fdml;
         totfdvl=totfdvl~fdvl;
	 totthmse=totthmse|lhltmse;
         endif;

/* run ols */

{ vnam,m,bols,stb,vcols,stderr,sigma,cx,rsq,resid,dbw } = ols(0,y,x);

/* calculate mse */

      yhat=x*bols;
      olsmse=sumc((y-yhat)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsum=bols;
         selsum=stderr;
      endif;
      if trial>=2;
         blsum=blsum~bols;
         selsum=selsum~stderr;
      endif;

/* obtain the first differences */

olsl=xl'*bols;
olsh=xh'*bols;
fdols=olsh-olsl;

      if trial==1;
         fdolsum=fdols;
	 totomse=olsmse;
      endif;
      if trial>=2;
         fdolsum=fdolsum~fdols;
	 totomse=totomse|olsmse;
      endif;

/* run ols (with logit trans) */

yl=ln(y./(1-y));

{ vnam,m,bolsl,stb,vcols,stderrl,sigma,cx,rsq,resid,dbw } = ols(0,yl,x);

/* calculate mse */

      yhatlt=exp(x*bolsl)./(1+exp(x*bolsl));
      ltmse=sumc((y-yhatlt)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsumo=bolsl;
         selsumo=stderrl;
      endif;
      if trial>=2;
         blsumo=blsumo~bolsl;
         selsumo=selsumo~stderrl;
      endif;

/* obtain the first differences */

olslo=exp(xl'*bolsl)./(1+exp(xl'*bolsl));
olsho=exp(xh'*bolsl)./(1+exp(xh'*bolsl));
fdolsl=olsho-olslo;

      if trial==1;
         fdolsuml=fdolsl;
	 tototmse=ltmse;
      endif;
      if trial>=2;
         fdolsuml=fdolsuml~fdolsl;
	 tototmse=tototmse|ltmse;
      endif;

/* run ols (with log trans) */

ly=ln(y+.01);

{ vnam,m,bolsll,stb,vcols,stderrll,sigma,cx,rsq,resid,dbw } = ols(0,ly,x);

/* calculate mse */

      yhatll=exp(x*bolsll)-0.01;
      ltmsel=sumc((y-yhatll)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsumol=bolsll;
         selsumol=stderrll;
      endif;
      if trial>=2;
         blsumol=blsumol~bolsll;
         selsumol=selsumol~stderrll;
      endif;

/* obtain the first differences */

olslol=exp(xl'*bolsll)-0.01;
olshol=exp(xh'*bolsll)-0.01;
fdolsll=olshol-olslol;

      if trial==1;
         fdsumll=fdolsll;
	 totmsell=ltmsel;
      endif;
      if trial>=2;
         fdsumll=fdsumll~fdolsll;
	 totmsell=totmsell|ltmsel;
      endif;

/* reset fails and add a trial for a successful trial */
      fails=0;
      trial=trial+1; 
      clear y;

   endif;

endo;

/* check that parameters are accurate */
print "descriptive stats for standard beta parameters";
meanc(bsum')~stdc(bsum'); 
print " ";
/* check the mse for the std beta */
  
print "standard beta mse=";;meanc(totbmse);
print " ";
print "alternative beta mse=";;meanc(totbamse);
print " ";
print "ols mse=";;meanc(totomse);
print " ";
print "logit trans ols mse=";;meanc(tototmse);
print " ";
print "het ols mse=";;meanc(totlhmse);
print " ";
print "logit trans het ols mse=";;meanc(totthmse);
print " ";
print "log trans het ols mse=";;meanc(totmsell);
print " ";

/* check mean and std. dev of first differences */

print"first differences for mean and variance for std beta (sim)";
meanc(totm')~stdc(totm');
meanc(totv')~stdc(totv');
print " ";

print "first differences for mean and variance for std beta (nosim)";
meanc(totmns')~stdc(totmns');
meanc(totvns')~stdc(totvns');
print " ";

print "first differences for mean and variance for alt beta (sim)";
meanc(totma')~stdc(totma');
meanc(totva')~stdc(totva');
print " ";

print "first differences for mean and variance for alt beta (nosim)";
meanc(totmnsa')~stdc(totmnsa');
meanc(totvnsa')~stdc(totvnsa');
print " ";

/* check for mean and std. dev. of first difference for OLS */
print"first differences for mean ols";
meanc(fdolsum')~stdc(fdolsum');
print " ";

/* check for mean and std. dev. of first difference for OLS (log trans) */
print"first differences for mean ols (logit trans)";
meanc(fdolsuml')~stdc(fdolsuml');
print " ";

print "first differences for mean and variance for het normal";
meanc(totfdmh')~stdc(totfdmh');
meanc(totfdvh')~stdc(totfdvh');
print " ";

print "first differences for mean and variance for logit trans.";
meanc(totfdml')~stdc(totfdml');
meanc(totfdvl')~stdc(totfdvl');
print " ";

/* check for mean and std. dev. of first difference for OLS (log trans) */
print"first differences for mean ols (log trans)";
meanc(fdsumll')~stdc(fdsumll');
print " ";

/* efficiency of expected values */
/* compare efficiency of standard beta with alternative beta */

effnba=((totmnsa-fdmt)^2);
effdba=((totmns-fdmt)^2);
effalt=100*(sqrt(sumc(effnba')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. alternative beta=";;effalt;
print " ";

/* compare efficiency of standard beta with ols */

effno=((fdolsum-fdmt)^2);
effdb=((totmns-fdmt)^2);
effols=100*(sqrt(sumc(effno')))./(sqrt(sumc(effdb')));
print "efficiency of beta vs. ols=";;effols;
print " ";

/* compare efficiency of standard beta with ols */

effnl=((fdolsuml-fdmt)^2);
effdb=((totmns-fdmt)^2);
effolsl=100*(sqrt(sumc(effnl')))./(sqrt(sumc(effdb')));
print "efficiency of beta vs. ols (log trans)=";;effolsl;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

effnbh=((totfdmh-fdmt)^2);
effdba=((totmns-fdmt)^2);
effhet=100*(sqrt(sumc(effnbh')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. het normal=";;effhet;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

effnbll=((totfdml-fdmt)^2);
effdbal=((totmns-fdmt)^2);
efflog=100*(sqrt(sumc(effnbll')))./(sqrt(sumc(effdbal')));
print "efficiency of standard beta vs. het normal (log trans)=";;efflog;
print " ";

/* compare efficiency of standard beta with log trans normal */

effnbl=((fdsumll-fdmt)^2);
effdba=((totmns-fdmt)^2);
efflogl=100*(sqrt(sumc(effnbl')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. OLS (log trans)=";;efflogl;
print " ";

/* efficiency for variance term */
/* compare efficiency of standard beta with alternative beta */

veffnba=((totvnsa-fdvt)^2);
veffdba=((totvns-fdvt)^2);
veffalt=100*(sqrt(sumc(veffnba')))./(sqrt(sumc(veffdba')));
print "efficiency of standard beta vs. alternative beta=";;veffalt;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

veffnbh=((totfdvh-fdvt)^2);
veffdba=((totvns-fdvt)^2);
veffhet=100*(sqrt(sumc(veffnbh')))./(sqrt(sumc(veffdba')));
print "efficiency of standard beta vs. het normal=";;veffhet;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

veffnbll=((totfdvl-fdvt)^2);
veffdbal=((totvns-fdvt)^2);
vefflog=100*(sqrt(sumc(veffnbll')))./(sqrt(sumc(veffdbal')));
print "efficiency of standard beta vs. het normal (log trans)=";;vefflog;
print " ";

/* check standard beta parameters' overconfidence */

npoca=(bsum-meanc(bsum'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesum^2;
dpoc=sqrt(sumc(dpoca'));
ocp=100*(npoc./dpoc);
print "standard beta parms' overconfidence=";;ocp;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check alternative beta parameters' overconfidence */

npoca=(bsuma-meanc(bsuma'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesuma^2;
dpoc=sqrt(sumc(dpoca'));
ocpa=100*(npoc./dpoc);
print "alternative beta parms' overconfidence=";;ocpa;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols parameters' overconfidence */

npoca=(blsum-meanc(blsum'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsum^2;
dpoc=sqrt(sumc(dpoca'));
oco=100*(npoc./dpoc);
print "ols parms' overconfidence=";;oco;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols (log. trans) parameters' overconfidence */

npoca=(blsumo-meanc(blsumo'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsumo^2;
dpoc=sqrt(sumc(dpoca'));
ocol=100*(npoc./dpoc);
print "ols (log. trans) parms' overconfidence=";;ocol;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols (log trans) parameters' overconfidence */

npoca=(blsumol-meanc(blsumol'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsumol^2;
dpoc=sqrt(sumc(dpoca'));
ocoll=100*(npoc./dpoc);
print "ols (log trans) parms' overconfidence=";;ocoll;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check heteroskedastic normal parameters' overconfidence */

npoca=(bsumh-meanc(bsumh'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesumh^2;
dpoc=sqrt(sumc(dpoca'));
ocph=100*(npoc./dpoc);
print "heteroskedastic normal parms' overconfidence=";;ocph;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check heteroskedastic logit transform parameters' overconfidence */

npoca=(bsuml-meanc(bsuml'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesuml^2;
dpoc=sqrt(sumc(dpoca'));
ocph=100*(npoc./dpoc);
print "heteroskedastic log. trans. parms' overconfidence=";;ocph;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check overconfidence for the first differences (mean) */

print "standard beta's mean first difference overconfidence";
nmoca=(totm-meanc(totm'))^2;
nmoc=sqrt(sumc(nmoca'));
dmoca=totms^2;
dmoc=sqrt(sumc(dmoca'));
ocm=100*(nmoc./dmoc);
ocm;
print " ";

/* check overconfidence for the first differences (variance) */

print "standard beta's variance first difference overconfidence";
nvoca=(totv-meanc(totv'))^2;
nvoc=sqrt(sumc(nvoca'));
dvoca=totvs^2;
dvoc=sqrt(sumc(dvoca'));
ocv=100*(nvoc./dvoc);
ocv;
print " ";

/* check overconfidence for the first differences (mean) */

print "alternative beta's mean first difference overconfidence";
nmocaa=(totma-meanc(totma'))^2;
nmoca=sqrt(sumc(nmocaa'));
dmocaa=totmsa^2;
dmoca=sqrt(sumc(dmocaa'));
ocma=100*(nmoca./dmoca);
ocma;
print " ";

/* check overconfidence for the first differences (variance) */

print "alternative beta's variance first difference overconfidence";
nvocaa=(totva-meanc(totva'))^2;
nvoca=sqrt(sumc(nvocaa'));
dvocaa=totvsa^2;
dvoca=sqrt(sumc(dvocaa'));
ocva=100*(nvoca./dvoca);
ocva;
print " ";

/* check percentage of parameters within 95% conf int. */

print "check standard beta's percentage within 95% conf int";
cib=abs(bsum-meanc(bsum')).<(1.96*sesum);
sumc(cib')/(trial-1);
print " ";

print "check alternative beta's percentage within 95% conf int";
ciba=abs(bsuma-meanc(bsuma')).<(1.96*sesuma);
sumc(ciba')/(trial-1);
print " ";

print "check ols percentage within 95% conf int";
cio=abs(blsum-meanc(blsum')).<(1.96*selsum);
sumc(cio')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cit=abs(blsumo-meanc(blsumo')).<(1.96*selsumo);
sumc(cit')/(trial-1);
print " ";

print "check het ols percentage within 95% conf int";
cih=abs(bsumh-meanc(bsumh')).<(1.96*sesumh);
sumc(cih')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cil=abs(bsuml-meanc(bsuml')).<(1.96*sesuml);
sumc(cil')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cill=abs(blsumol-meanc(blsumol')).<(1.96*selsumol);
sumc(cill')/(trial-1);
print " ";

/* get information about percentage of significant parameters */
/* mean values */
print "check percentage of sig parms (standard beta)";
tstatbm=(abs(totm./totms)).>1.96;
tstatbv=(abs(totv./totvs)).>1.96;
sumc(tstatbm')/(trial-1);
sumc(tstatbv')/(trial-1);
print " ";

print "check percentage of sig parms (alternate beta)";
tstatba=(abs(bsuma./sesuma)).>1.96;
sumc(tstatba')/(trial-1);
print " ";

print "check percentage of sig parms (basic ols)";
tstatbo=(abs(blsum./selsum)).>1.96;
sumc(tstatbo')/(trial-1);
print " ";

print "check percentage of sig parms (logit trans ols)";
tstatbl=(abs(blsumo./selsumo)).>1.96;
sumc(tstatbl')/(trial-1);
print " ";

print "check percentage of sig parms (het ols)";
tstatboh=(abs(bsumh./sesumh)).>1.96;
sumc(tstatboh')/(trial-1);
print " ";

print "check percentage of sig parms (het ols - logit trans)";
tstatblh=(abs(bsuml./sesuml)).>1.96;
sumc(tstatblh')/(trial-1);
print " ";

print "check percentage of sig parms (ols -- log tranform)";
tstatbm=(abs(blsumol./selsumol)).>1.96;
sumc(tstatbm')/(trial-1);
print " ";


/* obtain information about y's */
print "descriptive stats for y's for successful runs";
meanc(ystats);
print "descriptive stats for y's for failed runs";
meanc(ystatf);

/* obtain information about x's */
print "descriptive stats for x's for successful runs";
meanc(xstats');
print "descriptive stats for x's for failed runs";
meanc(xstatf');

print "number of failed runs";
crash;
output off;